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Abstract 

We present extensive numerical investigations of the structural relaxation dy- 
namics of a realistic model of the amorphous high-temperature ceramic a-Si3B3N7, 
probing the mean square displacement (MSD) of the atoms, the bond survival prob- 
ability (BSP), the average energy, the specific heat, and the two-point energy aver- 
age. Combining the information from these different sources, we identify a transition 
temperature T c ~ 2000 K below which the system is no longer ergodic and physical 
quantities observed over a time t b s show a systematic parametric dependence on 
the waiting time t w , or age, elapsed after the quench. The aging dynamics 'stiffens' 
as the system becomes older, which is similar to the behavior of highly idealized 
models such as Ising spin-glasses and Lennard-Jones glasses. 
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1 Introduction 

Physical aging of polymers and other 'soft' materials was studied experimentally al- 
most three decades ago by Struik and co-workers |2] . The interest of a wider community 
for this area was then greatly stimulated by work of Lundgren et al. [3] on the magnetic 
properties of spin glasses, and aging has since remained a source of important questions 
and insights in non-equilibrium statistical physics jl]. 

The basic features of aging were mainly established during the eighties in studies of 
magnetic linear susceptibility and autocorrelation function of spin glasses EH IB|, 
but are shared by many other systems, including, to cite only a few examples, type II 
superconductors P, glasses jTUl EEU U21 Ed Ej, granular materials [T3], and soft con- 
densed matter ^B]- Following a rapid quench of a relevant parameter, e.g. when the 
temperature of a glass is rapidly lowered below a threshold value, aging materials are 
unable to re-equilibrate within the available observation time. What is observed is a 
slow change of physical properties which occurs at a decelerating pace, and in a fashion 
largely independent of most details of the microscopic interactions. More precisely, physi- 
cal observables systematically depend on the system age t w , i.e. the time elapsed from the 
initial quench, and the observation time £ & s . Based on the relative values of t w and t t, s , 
two dynamical regimes can be identified: For t b s <C t w , macroscopic physical averages 
are approximately constant, and two-time correlations are related to the linear response 
through the fluctuation-dissipation theorem. Such a 'quasi-equilibrium' regime strongly 
resembles the time translational invariant dynamics of thermal equilibrium, except for the 
fact that many quantities, e.g. relaxation times, can carry a parametric age dependence. 
For t obs ^> t w , physical averages visibly drift and two-time averages no longer obey the 
fluctuation dissipation theorem. 

Recent improvements in time resolved spectroscopic techniques [d HH1 CHI ED] , have 
highlighted that, far from being a smooth flow, the drift of aging systems results from 
a highly irregular, so-called intermittent, process, whereby sudden and large configura- 
tional re-arrangements |2I| punctuate the prevailing equilibrium-like fluctuations. These 
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intermittent events may correspond to irreversible jumps from one metastable state to 
the next [221123] . 

Aging has raised a wealth of theoretical questions, for example, whether the non- 
equilibrium aging regime can be described through the concept of effective tempera- 
ture [2U 123 • Of direct relevance for the present investigations is the origin of the so- 
called pure aging behavior, i.e. the t Q b s /t w scaling of various observables in the aging 
regime [211 EEl and of modifyied scalings [2I| which are called sub-aging and super- or 
even hyper-aging. 

This work discusses numerical investigations of the dynamical and thermodynamical 
properties of a realistic model of a-SiaF^Ny for a large temperature range. The material is 
a representative of a class of nitridic ceramics with chemical composition a-SLB^Ni C|281 

EHl EU E21 EH EE E3 which can be synthesized via the sol-gel process |33|- 
These ceramics are of great technological interest and are presently considered for high- 
temperature engine applications. 

Nitridic ceramics appear to form amorphous covalent networks with a homogeneous 
distribution of the cations at least down to a scale of 1 nm[38j. Experimental studies 
have mainly focused on the quaternary compound a-SiBNsC. Under standard conditions, 
both a-SiBN 3 C and a-Si3B 3 N 7 are thermally stable and amorphous up to very high tem- 
peratures, e.g. about 2100 K for a-SiBNsC. They possess excellent mechanical and elastic 
properties: e.g. a-SiBN 3 C has a high bulk modulus of ca. 200 — 300 GPa, and is stable 
against oxidation up to 1700 K. 

The amorphous structure suggests a similarity with glassy materials, and hence the 
presence of aging properties of potential technological relevance. Probing the glass tran- 
sition of e.g. a-Si3B 3 N7, the basic representative of this class of nitridic ceramics, is im- 
possible at standard conditions, since decomposition takes place at T sa 1900 K, i.e. well 
before the ceramic melts. However, the melting transition might be accessible, as we have 
argued, at very high pressures [3Hj . 

The MC simulations presented below aim at the determination of the temperature 
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range below which the non-ergodic behavior of a-Si3B 3 N7 sets in, and at the characteri- 
zation of the aging dynamics in that low temperature regime. Admittedly, determining 
"for all practical purposes" whether a system is ergodic or not is not straightforward and 
is possibly an ill-posed questioned). As will become clear from the discussion, we refer 
to non-ergodic behavior as a collection of dynamical properties usually associated with 
glassiness, which, taken together, describe a strong deviation from thermal equilibrium. 

To detect the onset of ergodicity breaking and characterize the behavior in both the 
ergodic and non-ergodic regions, we consider a number of relevant dynamical and ther- 
modynamical measures: 1) The physical movement of the particles, as statistically char- 
acterized by the Mean Square Displacement (MSD), is diffusive in the ergodic region, and 
sub-diffusive elsewhere. 2) The related Bond Survival Probability (BSP)jH] drops to a 
finite average close to zero on a characteristic time scale at sufficiently high temperatures 
but becomes scale-free and age dependent in the non-ergodic region. There, it decreases 
as a logarithmic function of £ & s /i^, where a is close to one and depends on the type 
of particles forming the bond. 3) The specific heat Cy is calculated in three different 
ways, which agree in equilibrium but markedly differ in the non-ergodic region. 4) The 
two-time energy-energy average (fi(t w ,t b s ;T) defined in equation |H1 By construction, 
equals one in equilibrium or in a quasi-equilibrium regime, and is larger than one when the 
system is out of equilibrium. 5) The age dependence of the energies of the local minima 
or so-called inherent structures (IS) 02] characterizing the trajectories in the non-ergodic 
regime. These energies decrease logarithmically with the age. 

The rest of the paper is structured as follows: in Section El we introduce the model and 
simulation technique. In Section El we define the statistical measures used to probe the 
ergodicity breaking and present the corresponding results. Finally, Section 0] contains a 
broader discussion of the results, where we show that our findings are in general agreement 
with previous studies of the aging dynamics of much simpler atomistic models ^21 SHI IS] 
and also with experimental results in e.g. spin-glasses jlSj, polymer systems j2] and charge- 
density waves 
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2 Model and Techniques 

The model of a-SisBsNy consisted of 162 Si-atoms, 162 B-atoms and 378 N-atoms, 
respectively, in a 19.1 x 19.1 x 19.1 A 3 cubic box. As an interaction potential, we employed 
a two-body potential from the literature |17j based on fits to ab-initio energy calculations 
of crystalline and molecular compounds containing Si-, B- and N-atoms. This potential 
successfully reproduces experimental data regarding structural properties such as bond 
lengths, and bulk properties, e.g. vibrational frequencies, of the binary compounds SisN 4 
and BN and of molecules containing Si— N— B units. 

The starting configurations for our simulations were generated by relaxation from high- 
temperature melts. The simulations were performed at fixed temperature and volume, 
with a Monte-Carlo algorithm using the Metropolis acceptance criterion. In each update, 
an atom is randomly selected for an attempted move with a random direction, and with 
an average size chosen to achieve an acceptance rate of ~ 50 %. One Monte-Carlo 
cycle (MCC) corresponds to N atom = 702 such individual moves. Note that the kinetic 
energy (3/2/csT per atom) does not appear in MC-simulations, and that all quantities 
studied relate to the configurational energy. For comparison, we also performed MD 
simulations in the melt. Comparing mean square displacements determined from MD 
and MC simulations, respectively, we estimate that (at least in the high-temperature 
regime) one MCC corresponds to about 0.5 fs. At low temperature, the value of the 
microscopic time scale is not known precisely. However, this value is not crucial for the 
study of aging quantities which, as in our case, have a logarithmic time dependence or 
depend on the ratio of two times[48j. 

The temperatures investigated ranged from 25 to 10000 K; however, for T < 250 K 
no discernible dynamical evolution took place on the time scale of the simulations. For 
each temperature up to 3000 K and above 3000 K, 9 and 3 runs, respectively, of length 
ttotai = 2 x 10 5 MCC were performed. 

The energy E(t; T) as function of time was registered every 10 MCC. Along the in- 
dividual trajectories for T = 250, 7000 K, halting points xh were chosen, from which 
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conjugate gradient minimizations were performed. In the following, t init ~ 1000 MCC is 
the initialization time of the MC-simulations needed for the system to reach equilibrium 
in the ergodic regime (i.e. at high temperatures), while t w > ti n u is the waiting time before 
the observations begin. Furthermore, an ensemble of 100 runs of length t to tai = 10 6 MCC 
was studied, in order to investigate aging in greater detail, and with better statistics. 
Each starting configuration of the ensemble was generated by quenching the system from 
4000 K to 1250 K. 1 

The diffusivity of the atoms in the high temperature regime was calculated straight- 
forwardly from the mean square displacement as a function of time via the Einstein- 
Smochulowski relation (D oc M5D ft° 6s h . Finally, in the bond survival probability data, a 
B-N or Si-N bond was considered broken if the interatomic distance was above 2.0 A. 

In computing time averages and correlations we considered the observation interval 
/ = [t w ,t w + t b s ]. Time averages over I are denoted by an over-bar, while an angular 
bracket stands for an average over an ensemble of independent trajectories. 

3 Ergodicity-breaking and aging behavior 
3.1 Ergodicity breaking 

Evidence for ergodicity breaking can be extracted from a change of behavior in dy- 
namical quantities i.e. the MSD and BSP, as well as from thermodynamical quantities, 
i.e. the heat capacity, the energy average and the two-point energy average. All these 
data concur that ergodicity breaking occurs but differ to some degree with respect to the 
precise value of the transition temperature T c . The latter can be approximately located 
within the range 2000 - 2500 K. 

We consider the high temperature behavior first. For T > T' = 2100 K, the MSD 
between pairs of time-displaced atom positions 

MSD(t obs , t w ) = J2 [n(t w + Us) - n(t w )} 2 . (l) 

i 

1 Due to the rather expensive energy evaluation, performing the runs in this ensemble required ap- 
proximately a computing time of one year on a AMD/1800 PC. 
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only depends on the length of the observation interval t a b s = (t obs +t w ) —t w . The diffusion 
coefficients for the different atoms all have an Arrhenius temperature dependence, with 
activation barriers AE = 1.24 ± 0.09 eV, 1.27 ± 0.06 eV for silicon and boron atoms, 
respectively. 

In the same regime, the BSP is age independent and has a finite relaxation time tbsp- 
The values of tbsp were obtained by fitting the time dependence of the bond survival 
probability after the initialization time ti n n to a Kohlrausch- Williams- Watts law: [1H] 



yielding exponents (3 ~ 0.3 — 0.5. The time scale tbsp{T) has an Arrhenius-like tem- 
perature dependence, with activation barriers AE = 1.64 ± 0.02 eV and 1.60 ± 0.02 eV 
for Si-N and B-N bonds, respectively. Arrhenius plots for the diffusion coefficients of 
silicon and boron atoms and for the relaxation times tbs p are plotted in Fig. ^ A devi- 
ation from the Arrhenius scaling is apparent in both quantities as T approaches 2000 K 
from above. The activation barriers found for both diffusion and bond breaking are in a 
similar range, suggesting that a single Arrhenius time scale might characterize the high 
temperature dynamics. However, although diffusive processes need to be accompanied 
by bond-breaking, it is not clear whether this effect or steric blocking by neighbor atoms 
dominates the diffusion at various temperature scales. 

Deviations from activated behavior appear in the relaxation time as the temperature 
decreases. Approaching T' ~ 2100 K from above, the diffusion coefficients for B-, Si- and 
N-atoms all vanish for T — > T' as D oc (T — T') 7 with 7 = 1.6. A similar deviation of tbsp 
from an Arrhenius law is observed when T ~ 1720 K for B-N and T" « 1820 K for Si-N, 
respectively, are approached from above. Close to these temperatures, t B sp(T) grows 
roughly like a power law: tbsp(T) oc {T/T — l) s with 5b-n ~ —2.5, and Ssi-N ~ —2.3. 
At even lower temperatures, the particle motion becomes sub-diffusive and the BSPs 
fall off logarithmically. Both quantities become scale free but gain a characteristic age 
dependence, discussed in more detail in section 13.21 

Thermodynamical variables such as the specific heat Cy, the energy average, and a 
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Figure 1. Temperature dependence of relaxation times t B sp for Si-N (empty squares) 
and B-N (filled squares), and temperature dependence of diffusion coefficient D(Si) 
(empty circles) and D(B) (filled circles). Straight lines show fits to an Arrhenius law. 
The dotted vertical line indicates the critical temperature T c ~ 2000 K. 
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suitably normalized two-time energy-energy average require averaging over an ensemble. 
For a fixed amount of numerical effort, they are hence less precisely determined than 'local' 
dynamical properties, but offer nevertheless concurring insights into ergodicity breaking 
and aging behavior. 

The heat capacity was evaluated with three different prescriptions, all agreeing in 
thermal equilibrium. Note, that all specific heat data are reported on a per-atom basis. 
The disagreement observed at lower temperature shows that the system is non-ergodic for 
the observational time scale used. We note that agreement at higher temperatures does 
not neccessarily imply that the equilibrium state is reached jHO] > since such an agreement 
is only a neccessary but not a sufficient criterion for the system to be in equilibrium. 

Our first prescription is naively carried over from equilibrium thermodynamics. Leav- 
ing the T dependence of the mean energy understood on the r.h.s. of the equation, we 
define 

Cy(T) = (3) 

where the temperature derivative is performed numerically after the ensemble and time 
averaging denoted by angular brackets (()) and an overbar, respectively. 

Secondly, we consider simulations where the temperature is increased or decreased by 
a small amount ±AT pa ±0.1T at age t w . The trajectories are then time-averaged over 
the time t b s <C t w spent at T ± AT. 

This procedure yields 



rb m _ E(T + AT;t w )-E(T-AT;t w ) 

Finally, we gauge the variance of the energy fluctuations within the observation interval 
I — [t w ,t w + tabs] and estimate the heat capacity as 

= kBT 2 > (5) 

for a range of observation times £ & s which straddles t w . Figure where all Cy values are 
on a per-atom basis, shows that all three expressions agree for temperatures above 3000 
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K, but clearly disagree below 2000 K, which leads to a provisional value for the ergodicity 
breaking temperature T c « 2500 K. 

Within the non-ergodic region T < T c , we have studied C V (T) for t a b s <C t w , the 
quasi-equilibrium regime, and for t Q b s ^> t w the non-equilibrium regime. The dependence 
of Cy(T) on T, t b s and t w can also be obtained from the fitted time and age dependence 
of the mean energy given in Eq. below, by performing the time average and then 
differentiating with respect to T. For T < T c and t Q b s 3> t w , we find that C v has a 
small dependence on the length of the time interval t Q b s over which the time average is 
performed. This is connected to the logarithmic drift of the energy, which moreover is 
more pronounced just below T c than near T = (c.f. fig. EJ. 

C v mimics an experiment performed after some relatively long equilibration time 
t w ^> t b s and thus yields the most "realistic" value for the specific heat for all temper- 
atures. C v provides a way to check the validity of the fluctuation-dissipation theorem: 
it should coincide with the two other expressions above T c , and additionally, below T c in 
the restricted quasi-equilibrium region t Q b s < t w . Indeed, when increasing t a b s past t w the 
observed dynamics in C v changes from quasi-equilibrium to off-equilibrium in the tem- 
perature range between 2000 K and 3000 K (c.f. inset in fig. EJ)- Recall that the (trivial) 
contribution to Cy due to the change in kinetic energy is never included, i.e., in equilib- 
rium Cy — > 1.5&B for T — > K. Furthermore, the interaction potential V(r) vanishes for 
r — > oo, and thus this contribution to Cy goes to zero as T — ► oo. 

As shown in fig. the above prescriptions for the specific heat yield, as expected, 
almost identical results in the high-T ergodic dynamical regime, but differ at low T. 
This indicates that ergodicity is broken below T c w 2000 — 3000 K. Even in the low 
temperature region, C v and C v nevertheless roughly coincide in a restricted time interval 
which corresponds to the quasi-equilibrium regime t Q b s < t w . However, the inset clearly 
shows the large difference between C v in the quasi-equilibrium and the non-equilibrium 
regime, for temperatures below Tc- 
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Figure 2. Temperature dependence of the specific heats per atom C v (□), C v (o) and 
C v (•), for T G [500, 10000] K. For C v and C v , the waiting time f w and the observation 
time t bs are both equal to 10 5 MCC. C v is obtained by averaging over 10 equidistant 
t w values in the interval 10 5 — 2 • 10 5 , for each of three independent trajectories. The 
observation time was kept at t obs = 5 ■ 10 3 <C t w . All procedures agree for T > 3000K, the 
region in which the system equilibrates according to this measure. Inset: For temperatures 
T G [1000, 8000] K and waiting time t w = 10 4 , Cy is shown with symbols o and • 
corresponding to observation times t Q b s = 10 4 and t obs = 10 5 respectively. The high 
values (•) at low temperatures seen for t obs ^> t w are a consequence of the energy drift 
present in this dynamical regime. For the shorter observation time (o), the values are 
in good agreement with the C v values shown in the main panel with the same symbol 
(note the difference in scale). The agreement expresses the equivalence between drift and 
fluctuation within the quasi-equilibrium low temperature regime of the aging system. 
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Yet another probe of ergodicity-breaking is the two-time energy average defined by 

, (E(t w )E(t w + t obs )) 

0( ^'^ ;T) - (E(t w )E(t w )} ' (6) 

Above T c , and for t obs < t w below T c , we expect the dynamics to be time translational 
invariant, and hence = 1. Considering that the energy is negative and decreases mono- 
tonically with age, we expect an increase of with observation time, once the dynamics 
gains access to states of considerably lower energy. 

Figure El depicts for three temperatures, T = 1250 K, T = 1750 K and T = 2500 K, 
and for three values of the system age, t w = 3 ■ 10 3 , 10 4 and 10 5 MCC. The comparison 
clearly shows that ergodicity is broken as described at the lowest temperature (T = 1250 
K): due to the decrease of the configurational energy during aging, the two-time energy 
average increases as the observation time increases. At the intermediate temperature (T 
= 1750 K) the aging effect is still visible but difficult to quantify, and at the highest 
temperature (T = 2500 K) it has vanished completely. 

Taking all theses indicators, MSD, BSP, Cy b,c and into account, we conclude that 
a-SisB3N7 exhibits ergodicity breaking at T c « 2000 K. 

3.2 Aging properties 

In order to describe the non-ergodic regime in more detail we have studied the age 
dependence of: (i) the energy of the inherent structures and of the mean potential energy, 
(m) the bond survival probability (BSP), (m) the mean square displacement (MSD) of 
the particles, and (iv) the two point correlation function and the closely related energy- 
energy autocorrelation function. 

The ensemble-averaged energy (Eis{t w ,t)) of the inherent structures lying 'below' the 
trajectories observed in the interval [t w ,t w + t b s ] decreases logarithmically with t Q b s , as 
does the mean potential energy (E(t w ,t obs )). In both cases the logarithmic slope has a 
clear temperature dependence, but no significant t w dependence for t w > t init = 10 3 MCC. 
The energy simply decreases logarithmically with the simulation time elapsed after the 
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Figure 3. Energy self correlation functions 4>{t ohs ) for three different temperatures 1250 
K, 1750 K and 2500 K (bottom to top). Red, green and blue curves were recorded after 
waiting times t w — 3 • 10 3 , 10 4 and 10 5 MCC, respectively. The ensemble averages were 
calculated over 5 (1750 K), 6 (2500 K), and 100 (1250 K) different trajectories. Number 
inside the lower figure indicate waiting times t w (in MCC) after which measurements of 
4>(t obs ) began. 
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Figure 4. Time dependence of the average energies (E(t;T,x m i n )) = (Ejs(t,t w = 
Unit] T)) of the minima x min for selected temperatures T = 750 K, 1250 K, 1500 K, 1750 
K as function of time t = t init + t obs . 

quench, i.e. to a good approximation we find the functional form 

E(t obs ; T) = E (T) - A(T) In . ( 7 ) 

Here, E (T) = (E{t ini t;T)) and E (T) = (Eis(ti n i t ;T)) for potential energies E and 
minimum energies Ei$, respectively, and 

A IS (T) w 6.60 ■ 10~ 6 • T, A E w 6.64 ■ 10" 6 ■ T, (8) 

where the coefficients are measured in eV/(K-atom). Note that energies are reported on 
a per-atom basis. 

The logarithmic decay of the energies is shown in Fig. 01 and has previously been observed 
in other aging systems, see e.g. jTHl EJ E21 153| l5lj. 

Bond survival probabilities at T = 1250 K for Si-N and B-N bonds are plotted in 
Fig. El as a function of the scaled variable t/t® for three different ages t w = 10 4 , 10 5 , 8 • 10 5 
MCC. Even though the data collapse is not perfect, the scaling plots demonstrate that 
the bonds become, on average more resilient as the age increases. For the Si-N bond, the 

14 



1 




0.95 a o 



o 



0.9 



co 0.85 



CO 



0.8 



0.75 



0.7 



0.01 



0.1 



1 



10 



100 



1000 



obs' l w 



Figure 5. Aging dependence of the bond survival probability BSP for Si-N and B-N 
bonds, respectively, as function of t a b S /t°, at temperature T = 1250 K. 

best collapse is obtained for asiN = 1- Thus, the survival probability for the Si-N bond is 
a function of t/t w , a so-called pure aging behavior. The survival probability for the B-N 
bond, however, shows a sub-aging behavior with olbn — 0.75. For both types of bonds 
the survival probability decays in a logarithmic fashion after a nearly flat initial part, and 
no cut-off is reached at low temperatures, within the time scale of the simulations. 2 

Another quantity that exhibits an age dependence for T < Tq is the average mean 
square displacement of the atoms, given in Eq. ^ We observe a sub-diffusive behavior, 
MSD(t b s ,t w ) oc t obs with (3 < 1 {(3 = 1 corresponds to standard diffusion), in the quasi- 
equilibrium region (t obs <^t w ). 

The MSD remains sub-diffusive for t f, s > t w but with an increase in the exponent (3. 

In figure |6(a)| we show the mean square displacement MSD(t b s ,t w ) for three different 

waiting times t w , at temperature T = 1250 K. For t obs t w , we find (3 = 0.5, and the 
2 For temperatures different from 1250 K, the scaling exponents a can only be approximately deter- 
mined due to the small ensemble of simulation runs at these temperatures. The values for a lie in the 
range between 0.3 to 0.8 
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Figure 6. Double logarithmic plot of the mean square displacement as function of 
observation time t obs for different waiting times t w , at temperature T = 1250 K. The 
exponent (3 ~ 0.5 of the sub-diffusive motion MSD oc t ohs for waiting time t w is found 
by a fit of the data to a straight line, where only points up to the waiting time t w were 
included. (3 appears to be approximately independent of waiting time for t obs < t w . 
Although the collapse of the curves after t^-scaling in subfigure b) is not perfect, due to 
the fact that (3 starts increasing for t obs > t w , the data indicate strong sub-aging to be 
present. 

curves can be approximately collapsed (c.f. figure |6(b)| ) if t obs is scaled with The 
exponent a = 0.3 indicates a strong sub-aging behavior of the mean square displacement. 
This analysis was performed separately for all three types of atoms, Si, B and N. In all 
three cases, essentially the same results were found. 

The two-point correlation function (p(t w ,t obs ;T) (c.f. fig. EJ) remains very close to the 
equilibrium value 1 in the ergodic temperature region (T = 2500 K), and does not depend 
on waiting time. But in the glassy phase (T = 1750 K and T = 1250 K), 4>(t w ,t obs ]T) 
strongly deviates from its equilibrium value (0=1) for t obs > t w . Quasi-equilibrium 
behavior is only observed for observation times t obs < t w , that grow monotonically with 
t 3 

3 The closely related autocorrelation function CE{t w ,tobs',T) = (E(t w )-E(t w +t obs )) — (E(t w )) ■ (E(t w + 
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This monotonic dependence on t w of the time range during which quasi-equilibrium 
behavior is still observed, correlates with the stiffening of the response of the system 
characteristic for aging processes: The longer the system is allowed to equilibrate, the 
longer is the subsequent time range during which equilibrium-like behavior is observed. 
This effect concurs with the logarithmic decay of E(t; T) and (E(t; T)) (c.f. figure 0} an d 
fully agrees with the general aging behavior previously discussed. 

4 Discussion 

Numerical investigations of atomistic models with realistic interactions are computa- 
tionally very demanding, and for the same numerical effort, the accuracy of the results 
obtained is lower than what can be achieved for more simplified models. In spite of this 
limitation, we have shown how ergodicity is broken in a-Si3B3N7 and have characterized 
the dynamics in the low temperature phase, where the model behaves as a bona-fide aging 
system. In the aging regime the microscopic time scale of the dynamics is largely irrele- 
vant, and so is the precise value of the conversion factor between the time scales of Monte 
Carlo and molecular dynamics simulations. This assumption is supported by the fact that 
we have observed a similar slow drift in energy at low temperature MD simulations, i.e. 
in the non-ergodic regime, during studies of heat transport in a-Si3B 3 N 7 . jJS] Therefore, 
we trust that the qualitative and most of the quantitative aspects of our findings can be 
carried over to the real materials. 

The computational analysis of a-Si3B 3 N7 using MSD, BSP, C^ c and (f>(t w , t obs ; T) 
shows that this amorphous material can be expected to exhibit a glass transition with 
a concurrent break in the ergodicity at about T c « 2000 K if pressures high enough to 
prevent decomposition are applied. 4 

tobs)) also exhibits the expected aging behavior, i.e. Ce decreases to zero from an almost constant value 

(oc Cv(t w )) once t Q bs exceeds t w . Unfortunately, the data is too noisy for us to draw any quantitative 

conclusions, beyond the observation that the quasi-equilibrium regime grows monotonically with t w . 
4 Repeating these investigations for a large number of densities, we always find a similar freezing-in 

of the structure for T w T c . In other work, we have determined a critical point in the liquid-gas region 
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For T < T c , a-Si3B 3 N7 exhibits a general behavior in many respects similar to standard 
test systems for computer simulations (HUES El] (Lennard- Jones, a-Si0 2 ), insofar as we 
observe a freezing-in of the structure, a logarithmic decay of the energy, and a systematic 
waiting-time dependence of the bond survival probability, the mean square displacement 
of the particles and the two-time energy correlation function. 

The aging behavior of the heat capacity has been studied experimentally, for charge- 
density-wave systems (3ja], and for a spin-glass model [22] and discussed theoretically for 
two-level systems at very low temperatures [56 . 

The logarithmic drift towards lower energies occurring for T < T c is also observed in 
a number of model systems such as Lennard- Jones [HI] and soft-sphere glasses jS2]- This 
applies both to the actual trajectories and the time-sequence of observed local minima. 
For fixed simulation time, the deepest local minima are reached for T = 1750 K, which 
lies right below T c . Analogous observations are well-known from e.g. global optimization 
studies of complex systems, where it has been found that reaching the deepest local 
minima using Monte-Carlo-type search algorithms is achieved by spending most of the 
search time in the temperature interval slightly below the glass transition temperature |57| 
EH]- Thus this result provides another confirmation that ergodicity breaking takes place 
at T « T c . 

In summary, all the numerical evidence obtained from our simulations shows that the 
thermal quench of a S13B3N7 melt produces a material in a state of strong thermody- 
namic disequilibrium, which then relaxes through the aging process. As we have noted in 
the introduction, a-Si3B 3 N 7 is currently only synthesized via the sol-gel route. In earlier 
workjSH], we have shown that different synthesis routes (quench from melt, sol-gel, sinter- 
ing of nano-crystallites, etc.) result in amorphous compounds that exhibit a monotonic 

decrease in energy during simulations at T < 2000 K. This applies in particular to the 
of the ternary system (p cr w 1.3 GPa and, T cr w 8000K) [35]. Since the tendency to decompose is 
greatly reduced for a supercritical fluid, we predict that a-SiaBsNy should exhibit a glass transition at a 
temperature Tq ks 1700 — 2500 K and at a pressure of pc > 2 GPa. Up to now, high-pressure experiments 
have only been performed for T as 1000 K and p as 2.5 GPa. 
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relaxation of the structures that result from modelling the sol-gel synthesis route |oTH 16 lj. 
as we observed in both constant volume and constant pressure simulations [H2]- Thus, we 
surmise that aging phenomena analogous to those investigated here for a quenched melt 
also occur in the amorphous ceramic generated via the sol-gel route. 

The metastable states reached during aging become gradually lower in energy, which 
could lead to a strengthening of the material. On the other hand, the structural trans- 
formations which would accompany aging could lead to the formation of nano-crystallites 
and internal 'grain boundaries', whence the overall structure might develop e. g. cracks, 
and thus lead to a weakening of the mechanical properties of the material. This could 
especially be the case, if crystallites of binary or ternary phases are formed. However, 
for most envisaged applications, the time scales on which such processes might occur are 
beyond the accessible simulation time by many orders of magnitude. 
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